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Abstract 

I study the product of independent identically distributed D x D random 
probability matrices. Some exact asymptotic results are obtained. I find that 



J> . both the left and the right products approach exponentially to a probability 

00 . 

. matrix(asymptotic matrix) in which any two rows are the same. A parameter 

A is introduced for the exponential coefficient which can be used to describe 
ON . the convergent rate of the products. A depends on the distribution of indi- 



vidual random matrices. I find A = 3/2 for D = 2 when each element of indi- 



S 

. vidual random probability matrices is uniformly distributed in [0, 1] . In this 

case, each element of the asymptotic matrix follows a parabolic distribution 



function. The distribution function of the asymptotic matrix elements can be 

numerically shown to be non-universal. Numerical tests are carried out for a 

set of random probability matrices with a particular distribution function. I 

find that A increases monotonically from ~ 1.5 to oi 3 as D increases from 

3 to 99, and the distribution of random elements in the asymptotic products 

can be described by a Gaussian function with its mean to be i. 
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In recent years, there is an increasing interest in studying the properties of random 
matricesS. Random matrices can be used to describe disordered systems, chaos, and biolog- 
ical problems, and in the statistical description of complex nuclei^!. There are two different 
kind of problems related to random matrices. One of them is to study the statistical prop- 
erties of a single random matrix. WigneiH was the first to use a single large random matrix 
to explain the statistical behavior of levels in nuclear physics. The semicircular and the cir- 
cular theoremsB! had been obtained for the Gaussian unitary, Gaussian orthonormal, and 
Gaussian symplectic ensembles. The second problem is to study the statistical behavior of a 
product of random matrices. The product of random matrices has attracted much attention 
in most of recent works because many statistical problems in disordered systems and chaotic 
dynamical systems can be formulated as a study of such product. For the product of random 
matrices, there are beautiful Furstenberg§ and Oseleded theorems about the existence of 
Lyapunov characteristic exponents. However, a detailed analysis of the product of random 
matrices is very difficult due to the noncommutability of random matrices, and there are 
not so many general exact results about structures of the products of random matrices. 
Therefore, any exact results on the product of a particular type of random matrices should 
be interesting. 

In this work, I present some exact results of the product of independent identically 
distributed random probability matrices. A probability matrix T of D x D is defined as 
T(ij) > and ^2f =1 T(ij) = 1, i = 1,2, ...D, where T(ij) are the elements of T. Physically, 
such matrices can be used to describe the dynamics of a D state classical system if T t (ij) is 
interpreted as the probability for the system in i'th state jumping to j'th state at time t. If 
the hopping process is stochastic, the evolution of the system is described by the product of 
the random probability matrices T t . It is easy to show that the product of two probability 
matrices is still a probability matrix, ie. A = T{T 2 is a probability matrix provided T\ and 
T2 are two probability matrices of order D x D. In this paper, I show that the product of 
independent identically distributed random probability matrices approaches exponentially, 
in terms of the number of matrices in the product, to a matrix in which any two rows are the 



same. For 2x2 independent random probability matrices in which any element is uniformly 
distributed in [0, 1], I find that the asymptotic matrix elements have parabolic distribution. 
The main results are described by the following propositions. 

Proposition 1: Let {T k } be a set of D x D independent identically distributed ran- 
dom probability matrices in which all elements are random and have the same distribution 
function. If A(n) = YT k=1 T k = T n T n _i • • • T 2 Ti(the left product) and B(n) = Ul=iT k = 
T\T 2 • • • T„_iT n (the right product), then 



lim A(n) 



a(l) a(2) 
a(l) a(2) 



a(D) 
a(D) 



a(l) a(2) ... a(D) 



(1) 



and 



lim B(n) 



a(l) a(2) 
a(l) a(2) 



a(D) 
a(£>) 



(2) 



ya(l) a(2) ... a(D) y 

where a(i), (i = 1,2, . . . D), are positive random numbers with J2i a ij) — 1- However, the 
values of a's in the left product are fixed for a given sequence of random matrices while they 
keep changing in the right product. 

Before we prove this proposition, let us look at a special case of D = 2. Let 

f y n {\) l-y„(l)^ 
K y n {2) l-y n {2) J 



A{n) = T n • • • T 2 Tx 



(3) 



then, from A(n) = T n A(n — 1), we obtain the following recursion relations for y n (l) and 

J/n(2) 



y n {\) = x(l)y. 



n-l\ 



(l-a;(l))y n _i(2), 



(4a) 



y„(2) = x(2)y„_ 1 (l) + (1 - x(2)) 2/n _ 1 (2), 
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(4b) 



where x(l) and x{2) are the two independent random elements of T n , i.e., 



T 

-J- n. 



( x(l) l-x(l) ^ 
x(2) l-x(2) 



Therefore, 



which gives 



y n (l) - y n (2) = [x(l) - x(2)][y n _ 1 (l) - y n _ 1 (2)] 



S/n(l) - J/n(2) 



= |x(l) -x(2)| < 1. 



(5) 



S/n-i(l) - y„_i(2) 

Thus, we expect |y n (l) — 2/n(2)| approaches to zero exponentially. If x(l) and x(2) are 
uniformly distributed in [0, 1], then, for a large n, 



y n {\) - y n (2) ~ e" A ", 



(6) 



with 



A = - < Zn|x(l) - x(2)| >= - ln\x(l) - x(2)\dx(l)dx(2) = 3/2. 

Jo Jo 

Therefore, the left product of independent identically distributed 2x2 random probability 



matrices exponentially approaches to a matrix of the form 

a\ 1 — a\ 



\ 



(7) 



a 2 1 — a 2 

with A = 1.5. Similarly, it is easy to show that the same conclusion can be made for a right 
product. 

The above approach can be extended to the general cases. Without losing generality, 
we need only show that the values of elements in the first column of the product of random 
probability matrices approach to each other. Let 



A„ — 



Vnd) 



Vn{D) 



= T n A{n - 1). 



(8) 



We want to show that < \y n (k) — y n 0-)\ >— ex P , A > 0. It is not hard to see that 



D-l 



y n (i) ~ y n (D) = l x (^ k ) ~ x ( D i k )}[yn-i(k) - y n _i (£>)], 



(9) 



k=l 



where x(i,j) are matrix elements of T n which satisfies the conditions of a probability matrix. 
J2j x (h k) = 1 is used in the above derivation. Define 



Zn{i) = y n {i) - Vn{D) 



i = l,...,D-l, 



(10) 



equation (|9|) can be written in the following matrix form 



Zn{l) 
Zn{2) 



Z n {D 



C n 



J 



Zn-lW 



. z n -i(D - 1) i 



x(2,l) 



ar(D,l) 



a;(l,D-l)-a;(D,D-i; 
x(2,D- l)-x(D,D- 1] 



Zn-l(l) 
2n-l(2) 



x(D -1,1)- 1) • • • x(D - 1, D - 1) - x(L>, D-l] 



J V 



:n) 



where C n is a (£) — 1) x [D — 1) matrix related to probability matrix T n . It is not difficulty 
to show the following relation between T n and C n 



\T n - fjl\\ = (1 - /u)||C„ - 



(12) 



that is the eigenvalues of matrix C n are D — 1 of eigenvalues of matrix T n (one of the 
eigenvalues of T n with value 1 is excluded)!. It is well known that the magnitudes of 
eigenvalues /ij, i — 1, . . . , D, of a D x D probability matrix are not greater than 1, i.e. 
< 1. Furthermore, matrices C's do not have any common eigenvectors^. Therefore, 



< 



z n \l) 

Zn{2) 



Z n {l) 



-XN 



(for a large N) 



(13) 



with A > 0, i.e., the values of elements in the first column of the left product A approach 
exponentially to each other. It is easy to show that the same result is true for any other 
column of A. Thus, relation ([I]) holds. Similarly, it can be shown that relation (|2|) also 
holds. 

This result is not really surprising. It is known that the magnitudes of eigenvalues of a 

probability matrix are equal to or smaller than 1. A probability matrix always has an eigen- 

/1\ 



value 1 with the corresponding eigenvector(mode) 



(There may exist other eigenvectors 



with eigenvalue 1). The eigenvalues of the product will either approach to or stay at 1 
when such matrices are multiplied together. Because each matrix is random and indepen- 
dent, these matrices are not commutable among themselves, and they do not in general have 



the same eigenvector except 



which will keep unchanged since its eigenvalue is equal to 



\1/ 



1. Therefore, all other modes are mixed together, and decay with the multiplication. The 
relation ([!]) is then expected. 

Proposition 2: Let {T^} be a set of 2 x 2 independent identically distributed random 
probability matrices in which all elements are uniformly distributed in [0, 1], proposition 1 
guarantees that 



lim A(n) = lim T n T n _i • • • T 2 T X 



and 



lim Bin) = lim T{T 2 ■ ■ ■ T n ^{T n 



1 i \ 

a 1 — a 



a 1 — a 



a 1 — a 



a 1 — a 



Then a is a random variable whose distribution function is 



f(a) = 6a(l — a). 
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(14) 



To prove proposition 2, we notice that a is a random variable which obeies recursion 
relation 



a n = x(l)a n -i + x(2)(l-a n - 1 ), (15) 

where x(l) and x(2) are two indenpendent random numbers uniformly distributed in [0, 1]. 
Since f(a) is the asymptotic distribution function of a n , i.e., n — > oo, a n and a n _i should 
have the same distribution function f(a) when n — > oo. Therefore, we can obtain the 
following equation in the integral form for distribution function f(a) 

f( a ) = f(b)5(a — bxi + bx 2 — X2)dxidx 2 db. (16) 

Jo Jo Jo 

Substitute 5 (a — bx\ + bx 2 — x%) by 

/oo 
e iq { a-bx 1+ bx 2 -x 2 ) dq ^ 

and integrate over xi and x 2 , equation (|IED becomes 

f{a) = 2~Jo db b(b^T) L dq ? • (18) 



2^ 



b(b - 1) J-oo * q 2 

Differentiate equation flTif ) with respect to a twice and notice f(a) = /(l — a), we can show 
that f(a) satisfies differential equation 

f" = ~77^r/. ( 19 ) 

a(l — a) 

with boundary conditions /(0) = /(l) = 0. Equation (0) can be easily solved by power- 
series expansion method since a = (or a = 1) is a regular singular point. The solution of 
this equation (normalize to 1) is that of equation fll4|) . 

I have carried out some numerical simulations to further confirm the results in the above 
propositions. Figure [I] is the distribution function of an element of the asymptotic matrix 
of the product of 2 x 2 independent random probability matrices in which elements are 
uniformly distributed in [0,1]. Solid line is the numerical result, and dashed line is the 
analytical expression (|14[) . They agree very well with each other. In order to check whether 



the distribution function of the product depends on the distribution function of individual 
random probability matrices, I also study the product of independent random probability 
matrices in which elements are distributed in [0, 1] according to function 



Distribution (|20|) is chosen because it is easy to generate on a computer. Although I cannot 
find the distribution function for the matrix elements of the asymptotic product analytically 
in this case, numerical results can easily be obtained. Figure |] is the single variable distri- 
bution functions of random matrix elements in the right product of such random probability 
matrices of order of 2 x 2 and 4x4. The dashed lines are the numerical results, and the 
solid lines are the fits of Gaussian functions. The numerical results can be well described 
by a Gaussian function with its mean to be j^. Compare with that in figure p], we can see 
that the distribution function of the product of independent identically distributed random 
probability matrices depends on the distribution of individual random matrices. In other 
words, unlike the large number theorem for random numbers, the distribution function of 
the product is not universal. 

I check numerically the results in proposition 1 by using the random probability matrices 
whose elements are independently, except the constrains of a probability matrix, distributed 
in [0, 1] according to equation (f20D. I compute the decay of < 1) — A n (2, 1)| > with 

n, where < . . . > denotes ensemble average, and A n (i,j) are the elements of the product of 
n random matrices. Figure |3| is In < |A n (l, 1) — A n (2, 1)| > vs. Inn for the left product of 
3 x 3, 4 x 4, 8 x 8, 16 x 16, and 99 x 99 random matrices. 100 ensembles are used in the 
numerical study. Figure |] is a similar plot (as figure 1) for a right product. The exponential 
decay of the quantity is clearly shown in these figures. Numerically, I find that A increases 
monotonically from ~ 1.5 to ~ 3 as D increases from 3 to 99. 

In conclusion, I have shown that both left and right products of a sequence of independent 
identically distributed random probability matrices exponentially approach to a probability 
matrix in which all elements in any column vector are the same. An exponential exponent 




(20) 



S 



is used to describe this approach rate. I also find that A increases monotonically from ~ 1.5 
to ~ 3. as D increases from 3 to 99 when the distribution function of individual random 



matrices is described by equation pD). I also find that A = 3/2 for D = 2 when random 
matrix elements are uniformly distributed in [0, 1]. It is well known that at least one of the 
eigenvalues of a probability matrix is equal to 1 while the rest of them are distributed in 
[0, l]i. A large D means that the product has more channel to decay to the stable structure 
(H) or Thus it is expected that the decay is faster for large D, ie., A increases with D. In 
order to understand the meaning of the results, let us look at a physical model system of D 
states. Assume the system can move randomly from i'th state to j'th state with probability 
T t (ij) at time step t. If the system starts from an initial distribution, one might want to know 
the probability of the system in i' state, i.e., the distribution function, after a long time. The 
questions may be whether there is a stable distribution (equilibrium state), and/or what it is 
if there is one. Obviously, the long time distribution(s) are the non-trivial left eigenstate(s) 
of the right product of the random matrices T t . In equation (Q), although the structure for 
a product will not change when n is larger than certain value, elements a(i) do change as 
another independent random probability matrix is multiplyed to the product. Therefore, the 
system does not have a stable distribution as one expected since the dynamics of the system 
is a stochastic process, and transition probabilities keep change with time. Both of the 
right and left products, however, have a non-trivial unique right eigenstate with eigenvalue 
1. The eigenstate takes the same value in each of its D components. Unfortunately, I 
am not able to obtain any meaningful non-trivial results by applying the propositions to 
this simple model system. It will be interesting to find some interesting physical systems 
in which the propositions can be used to extract useful information. In contrast to the 
the sum of independent random variables whose distribution is Gaussian no matter what 
is the distribution function of individual random variables, the distribution function of a 
element in the product of independent identically distributed random matrices depends on 
the distribution of individual random matrices. Therefore, the distribution function of the 
product of independent identically distributed random matrices is not universal. 
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FIGURES 

FIG. 1. Distribution function of a random matrix elements in the right product of independent 
identically distributed 2x2 random probability matrices in which all elements are uniformly 
distributed in [0, 1]. The dash line is the numerical result, and the solid line is f(a) = 6a(l — a). 

FIG. 2. Distribution function of an arbitrary random matrix element in the right product of 
independent identically distributed 2 x 2, 4 x 4 random probability matrices in which all elements 
are distributed in [0, 1] according to equation (|20|). The dashed lines are the numerical results, and 
the solid lines are the fits of Gaussian functions. 

FIG. 3. In < |j4 n (l,l) — A n (2, 1)| > vs n of the left product of independent identically 
distributed random matrices of 3 x 3, 4 x 4, 8 x 8, 16 x 16, 99 x 99, with 100 ensembles. A's can be 
obtained from the slopes. The slopes increase monotonically from A ~ 1.5 to A ~ 3. as D changes 
from D = 3 to D = 99. 

FIG. 4. In < |^4 n (l,l) — A n (2, 1)| > vs n of the right product of independent identically 
distributed random matrices of 3 x 3, 4 x 4, 8 x 8, 16 x 16, 99 x 99. A's can be obtained from the 
slopes. The slopes increase monotonically from A ~ 1.5 to A ~ 3. as D changes from D = 3 to 
D = 99. 
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